clear; clc; close all;

%% 参数设置
k = 3;
h = 0.01;
T = 600;
transient = 550;
N = round(T / h);
N_trans = round(transient / h);

% 参数空间划分
mu_vals = linspace(-3*pi, 3*pi, 100);
eta_vals = linspace(-3*pi, 3*pi, 100);

% 结果矩阵
state_map = zeros(length(eta_vals), length(mu_vals));

max_period = 6; % 最大检测周期

%% 参数空间扫描（可用parfor加速）
parfor i_eta = 1:length(eta_vals)
    eta = eta_vals(i_eta);
    local_state_row = zeros(1, length(mu_vals));
    for i_mu = 1:length(mu_vals)
        mu = mu_vals(i_mu);
        x = [1e-9; 0; 0; mu; eta];
        f = @(x) system_ode(x, k);

        % 丢弃暂态
        for t = 1:N_trans
            x = RK4_step(f, x, h);
        end

        % 采集后期轨迹y变量
        y_seq = zeros(N - N_trans, 1);
        for t = 1:(N - N_trans)
            x = RK4_step(f, x, h);
            y_seq(t) = x(2); % y变量
        end

        % 分类
        local_state_row(i_mu) = classify_P_CH(y_seq, max_period);
    end
    state_map(i_eta, :) = local_state_row;
end

%% 绘图
figure;
imagesc(mu_vals/pi, eta_vals/pi, state_map);
set(gca,'YDir','normal');

% 颜色映射，8类
cmap = [
    1.0, 1.0, 0.0;  % P0 黄
    0.0, 0.0, 0.0;  % P1 黑
    0.0, 0.0, 1.0;  % P2 蓝
    0.0, 1.0, 1.0;  % P3 浅蓝
    0.0, 1.0, 0.0;  % P4 绿
    0.0, 0.5, 0.0;  % P5 深绿
    0.5, 0.0, 0.5;  % P6 紫
    1.0, 0.0, 0.0;  % CH 红
];
colormap(cmap);

colorbar('Ticks',0:7, ...
    'TickLabels',{'P0','P1','P2','P3','P4','P5','P6','CH'});
caxis([-0.5 7.5]); % 限制颜色条范围

xlabel('\mu / \pi');
ylabel('\eta / \pi');
title('Fig.5(a) 吸引子盆地图 - 周期与混沌分类');

xticks(-3:1:3);
yticks(-3:1:3);
grid on;

% 统计各类别数量
counts = zeros(1,8);
total_points = numel(state_map);

for c = 0:7
    counts(c+1) = sum(state_map(:) == c);
end

% 计算占比（百分比）
percentages = counts / total_points * 100;

% 类别标签
labels = {'P0','P1','P2','P3','P4','P5','P6','CH'};

fprintf('各状态占比：\n');
for i = 1:8
    fprintf('%s: %.2f%%\n', labels{i}, percentages(i));
end

% 各状态占比：
% P0: 3.60%
% P1: 44.80%
% P2: 11.28%
% P3: 2.52%
% P4: 2.28%
% P5: 2.24%
% P6: 0.48%
% CH: 32.80%

%% 子函数：系统方程
function dx = system_ode(x, k)
    dx = zeros(5,1);
    dx(1) = x(2) + x(3) - k * cos(x(5)) * x(2);
    dx(2) = -x(1) + x(3);
    dx(3) = -x(1) - x(3) + k * cos(x(4)) * x(1);
    dx(4) = x(1);
    dx(5) = x(2);
end

%% 子函数：RK4一步积分
function x_next = RK4_step(f, x, h)
    k1 = f(x);
    k2 = f(x + h/2 * k1);
    k3 = f(x + h/2 * k2);
    k4 = f(x + h * k3);
    x_next = x + h/6*(k1 + 2*k2 + 2*k3 + k4);
end

%% 子函数：周期+混沌分类函数
function class_num = classify_P_CH(traj, max_period)
    % 取峰值
    [pks, locs] = findpeaks(traj);

    if length(pks) < 2
        class_num = 0; % P0 定点
        return;
    end

    % 归一化峰值
    pks_norm = pks / max(abs(pks));

    for period = 1:max_period
        repeated = true;
        for i = 1:(length(pks_norm)-period)
            if abs(pks_norm(i) - pks_norm(i+period)) > 0.1
                repeated = false;
                break;
            end
        end
        if repeated
            class_num = period; % 周期判定
            return;
        end
    end

    % 无周期判定则为混沌
    class_num = 7; % CH
end
